## 在这个项目中主要存在4个可以直接优化点，分别是：
# 1、迭代次数 N_sim，如果修改 N_sim 则需要相应的改变图的位置。
# 2、fitness, 表达式或者系数都可以调整。
# 3、突变比例

getwd()

rm(list = ls())
seed <- 3
set.seed(seed)

# install.packages('ggplot2')

library(mgcv) 
library(ggplot2)
library(reshape2) 
library(tidyverse)
source("toolbox.R")


######################################
### settings for simulation number ###
######################################

# parameter
mech <- 5 # choose 1,2,3,4 for other mechanism
N_sim <- 5000  # 迭代次数

# a_list <- c(0,0.1,0.2,0.3,0.35,0.4,0.5,0.6,0.7,0.8,0.9)
# a_list <- c(0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42)
# alpha_list <- c(-0.5,-0.3,0,0.3,0.5)
# beta_list <- c(-0.5,-0.3,0,0.3,0.5)

# df <- expand_grid(a_list,alpha_list,beta_list)

# (nrow(df)/5*4+1):(nrow(df)/5*5))

for (i in 2:11) {
  # df <- as.data.frame(df)
  # a <- df[i,1]
  # b <- 1-a
  # alpha <- df[i,2]
  # beta <- df[i,3]

  a <- 0.3
  b <- 1-a
  alpha <- 0
  beta <- 0
  
  print(c(a,b,alpha,beta))
  
  # 描述不同fitness定义方式
  # for antagonist
  F_antagonist <- function(mech){
    E1 <- (1-H_V_A) # n=1 for A_min & P_max
    E2 <- (1-H_V_A)-alpha*(1-H_V_Po)
    E3 <- (1-H_V_A)
    E4 <- (1-H_V_A)+beta*(1-H_V_Po)
    E5 <- (1-H_V_A)+beta*(1-H_V_Po)
    
    E <- list(E1,E2,E3,E4,E5)
    return(E[[mech]])
  }
  
  # for pollinators
  F_pollinators <- function(mech){
    E1 <- (1-H_V_Po)   # n=1 for A_min & P_max
    E2 <- (1-H_V_Po)-alpha*(1-H_V_A)
    E3 <- (1-H_V_Po)-alpha*(1-H_V_A)
    E4 <- (1-H_V_Po)-alpha*(1-H_V_A)
    E5 <- (1-H_V_Po)+alpha*(1-H_V_A)
    
    E <- list(E1,E2,E3,E4,E5)
    return(E[[mech]])
  }
  
  # for plant
  F_plant <- function(mech){
    E1 <- (a*H_A_V+b*(1-H_Po_V))
    E2 <- (a*H_A_V+b*(1-H_Po_V))
    E3 <- (a*H_A_V+b*(1-H_Po_V))
    E4 <- (a*H_A_V+b*(1-H_Po_V))
    E5 <- (a*H_A_V+b*(1-H_Po_V))
    
    E <- list(E1,E2,E3,E4,E5)
    return(E[[mech]])
  }
  
  ###################################################
  ## load field observation matrix AP_obs, PV_obs ###
  ###################################################
  
  # # data 应为实际观测值
  # AP_obs  <- matrix(data = NA,nrow = 48,ncol = 11) # 这里设置植食昆虫有50种，榕树45种，有机化合物244种
  # PoP_obs <- matrix(data = NA,nrow = 18,ncol = 11)
  # PV_obs  <- matrix(data = NA,nrow = 11,ncol = 22)
  
  AP_obs  <- read.csv("data/AP_Africa_B.csv", header = TRUE, as.is = TRUE, row.names = 1)
  PoP_obs <- read.csv("data/PoP_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)
  PV_obs  <- read.csv("data/PV_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)
  
  
  ############################
  ### Preparation:Sample #####
  ############################
  
  index_choose <- c(35,40,16,40,74,12,14,15,1,19) # 按照原则一抽取的 n 种 VOCs
  PV_obs_1 <- PV_obs[,c(index_choose)]
  index_all <- 1:ncol(PV_obs)
  
  PV_obs_2 <- PV_obs[,c(index_all[-index_choose])]
  PV_obs_3 <- PV_obs_2[,c(colSums(PV_obs[,index_all[-index_choose]])<=i)]
  
  index_1 <- 1:ncol(PV_obs_3)
  # index_ <- sample(index_1 , 56)
  # print(index_)
  # PV_obs_4 <- PV_obs_3[,c(index_ )]
  
  PV_obs <- cbind(PV_obs_1,PV_obs_3)
  
  ########################################
  
  nA  <- nrow(AP_obs) # 植食动物的数量 （行计数）
  nP  <- ncol(AP_obs) # 植物的数量 （列计数）
  nV  <- ncol(PV_obs) #  VOC 的数量
  nPo <- nrow(PoP_obs) # 传粉蜂的数量
  
  ## 生成模拟的初始化矩阵
  # 1、初始矩阵不影响模拟结果
  # 2、初始时要保证任意植物至少有与一种植食动物有关；任意VOC都至少由一种植物释放。因此下述的colSum 的所以值都大于 0
  AP  <- matrix(rbinom(nA*nP,1,0.5),nA,nP)
  PV  <- matrix(rbinom(nP*nV,1,0.5),nP,nV)
  PoP <- matrix(rbinom(nPo*nP,1,0.5),nPo,nP)
  colSums(AP) #  确保 colSum 的所以值都大于 0
  colSums(PV) 
  colSums(PoP)
  # 每一次循环基因变异的比例，该值仅影响收敛速度，不影响结果。
  M1 <- 0.2*nP*nV/2 # PV 矩阵突变的比例
  M2 <- 0.2*nA*nP # AP 矩阵突变的比例
  M3 <- 0.2*nPo*nP # PoP 矩阵突变的比例
  
  #######################################
  ### Analyses1: Simulation process #####
  #######################################
  # create plot
  # par(mar=c(4,4,2,9))
  png(
    filename = paste0("grid_search_output/date_12_9_seed_",seed,"_",i,".png"), # 文件名称
    width = 650,           # 宽
    height = 480,          # 高
    units = "px",          # 单位
    bg = "white",          # 背景颜色
    res = 72)              # 分辨率
  
  par(mar=c(4,4,2,9))
  plot(-1, xlim = c(0,N_sim), ylim = c(0,1), ylab = "Fitness", xlab = "Time")
  legend(x=N_sim+N_sim/40, y=0.4, title = "Simulated", legend = c("Plant Fitness","Antagonist Fitness","pollinators Fitness"),
         pch = c(3,4,5), col = c("green3","red","blue"), xpd=TRUE, bty="n",title.font = 2) 
  
  # legend(x=N_sim+N_sim/40, y=0.8 , title = "Parameter", legend = c(paste0("a = ",a), paste0("b = ",b), paste0("alpha = ",alpha), paste0("beta = ",beta)),
  #        xpd=TRUE, bty="n",title.font = 2) 
  
  
  # 根据实地观察，计算并绘制熵值和fitness
  {
    AV_obs  <- as.matrix(AP_obs) %*% as.matrix(PV_obs)
    PoV_obs <- as.matrix(PoP_obs) %*% as.matrix(PV_obs)
    # "H_A_VOC" 函数 i:矩阵 o:互信息/条件熵，toolbox
    # H_V   <- H_A_VOC(PV_obs)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV_obs)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV_obs)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV_obs)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV_obs)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP_obs)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP_obs)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV_obs)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV_obs)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP_obs)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP_obs)[["Hn_V_S"]]
    
    # 对抗系统中的 fitness 表达, toolbox
    E_plant       <- F_plant(mech)
    E_antagonist  <- F_antagonist(mech)
    E_pollinators <- F_pollinators(mech)
    
    
    
    # E_obs <- c(0, E_plant, E_antagonist, E_pollinators, H_V,
    #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    # names(E_obs) <- c("N","E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
    #
    ## 野外实际观测值划线，图A
    # 条件熵 for PV, AV and AP
    # abline(h = H_P_V, lty = "dotted", pch = 8)
    # abline(h = H_A_V, lty = "solid", pch = 8)
    # abline(h = H_A_P, lty = "dashed", pch = 8)
    ## 植物和植食动物的fitness，图B
    abline(h = E_antagonist, col = "red")
    abline(h = E_pollinators, col = "blue")
    abline(h = E_plant, col = "green3") # the same as H(A|V) in mech = 1
  }
  
  
  # 模拟过程
  A_PV  <- list()  # 每次模拟后存储 PV矩阵 的 list
  A_AP  <- list()  # 每次模拟后存储 AP矩阵 的 list
  A_PoP <- list()  # 每次模拟后存储 PoP矩阵 的 list
  # E = matrix(NA, nrow = N_sim +1, ncol = 15) # E 用来存储每个模拟之后的变量值 N*14
  # colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
  E = matrix(NA, nrow = N_sim +1, ncol = 8) # E 用来存储每个模拟之后的变量值 N*14
  colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators",  "H_A_V", "H_V_A","H_Po_V","H_V_Po")
  
  
  ## 基于随机初始化的第一次模拟值
  n = 1
  {
    AV  <- AP %*% PV
    PoV <- PoP %*% PV
    # "H_A_VOC" 在toolbax.R 中输入矩阵，输出互信息和条件熵
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant       <- F_plant(mech) 
    E_antagonist  <- F_antagonist(mech) 
    E_pollinators <- F_pollinators(mech) 
    
    # 绘图
    # points(n,  H_P_V, col = "red", pch = 24) 
    # points(n,  H_A_V, col = "red", pch=  21)
    # points(n,  H_A_P, col = "red", pch=  22) 
    # 图B
    points(n,  E_pollinators, col = "blue", pch = 5)
    points(n,  E_antagonist, col = "red", pch = 4)
    points(n,  E_plant, col = "green3", pch = 3)
    
    #写入 E, A_AP, A_PV 的第N次
    # E[n,] <- c(n, E_plant, E_antagonist, E_pollinators, H_V,
    #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    E[n,] <- c(n, E_plant, E_antagonist, E_pollinators,
               H_A_V, H_V_A, H_Po_V, H_V_Po)
    
    A_AP[[n]] <- AP
    A_PV[[n]] <- PV
    A_PoP[[n]] <- PoP
  }
  
  
  ### 植食动物和植物依次使其 fitness 最大化 迭代模拟
  for (n in 1:(N_sim/4)) # N_sim : 模拟次数。由于模拟过程包括植食动物、传粉蜂的基因突变和植物两次基因突变四个部分，因此需要 "/4"
  {
    print(n)
    ## 传粉蜂：最大化fitness
    for (m in 1:M3) {
      PoP_new <- random.sample_1element1(PoP) #function details see toolbox.r
      PoV <- PoP_new %*% PV
      # H_V   <- H_A_VOC(PV)[["Hn_V"]]
      # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
      # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
      H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
      H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
      # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
      # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
      H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
      H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
      # H_Po_P <- H_A_VOC(PoP_new)[["Hn_S_V"]]
      # H_P_Po <- H_A_VOC(PoP_new)[["Hn_V_S"]]
      
      
      # fitness function
      E_plant_new       <- F_plant(mech) 
      E_antagonist_new  <- F_antagonist(mech) 
      E_pollinators_new <- F_pollinators(mech) 
      
      # 基于优化机制判断是否保留此次基因突变
      if(E_pollinators_new > E_pollinators)
      {E_antagonis  <- E_antagonist_new
      E_plant       <- E_plant_new
      E_pollinators <- E_pollinators_new
      PoP <- PoP_new 
      }
    }
    
    # 更新矩阵，并计算互信息量和条件熵
    PoV <- PoP %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant       <- F_plant(mech) 
    E_antagonist  <- F_antagonist(mech) 
    E_pollinators <- F_pollinators(mech) 
    
    # 写入模拟完植食动物的 E, A_AP, A_PV
    # E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators, H_V, 
    #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    
    E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators,
                   H_A_V, H_V_A, H_Po_V, H_V_Po)
    
    A_AP[[4*n-2]] <- AP
    A_PV[[4*n-2]] <- PV
    A_PoP[[4*n-2]] <- PoP
    
    # # 绘图
    # points((4*n),  H_P_V, col = "red", pch = 24) 
    # points((4*n),  H_A_V, col = "red", pch=  21)
    # points((4*n),  H_A_P, col = "red", pch=  22)
    # 图B
    points((4*n),  E_pollinators, col = "blue", pch = 5)
    points((4*n),  E_antagonist, col = "red", pch = 4)
    points((4*n),  E_plant, col = "green3", pch = 3)
    
    ## 被子植物：最大化fitness
    for (m in 1:M1) {
      PV_new <- random.sample_1element1(PV) #function details see toolbox.r
      PoV <- PoP %*% PV_new
      AV  <- AP %*% PV_new
      # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
      # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
      # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
      H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
      H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
      # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
      # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
      H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
      H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
      # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
      # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
      
      
      # fitness function
      E_plant_new       <- F_plant(mech) 
      E_antagonist_new  <- F_antagonist(mech) 
      E_pollinators_new <- F_pollinators(mech) 
      
      # 基于优化机制判断是否保留此次基因突变
      if(E_plant_new > E_plant)
      {E_antagonis  <- E_antagonist_new
      E_plant       <- E_plant_new
      E_pollinators <- E_pollinators_new
      PV <- PV_new
      }
      
    }
    
    # 更新矩阵，并计算互信息量和条件熵
    AV  <- AP %*% PV
    PoV <- PoP %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant       <- F_plant(mech) 
    E_antagonist  <- F_antagonist(mech) 
    E_pollinators <- F_pollinators(mech) 
    
    # # 写入模拟完植食动物的 E, A_AP, A_PV
    # E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators, H_V, 
    #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators,
                   H_A_V, H_V_A, H_Po_V, H_V_Po)
    A_AP[[4*n-1]] <- AP
    A_PV[[4*n-1]] <- PV
    A_PoP[[4*n-1]] <- PoP
    
    # # 绘图
    # points((4*n),  H_P_V, col = "red", pch = 24) 
    # points((4*n),  H_A_V, col = "red", pch=  21)
    # points((4*n),  H_A_P, col = "red", pch=  22)
    # 图B
    points((4*n),  E_pollinators, col = "blue", pch = 5)
    points((4*n),  E_antagonist, col = "red", pch = 4)
    points((4*n),  E_plant, col = "green3", pch = 3)
    
    
    ## 植食动物：最大化fitness
    for (m in 1:M2) {
      AP_new <- random.sample_1element1(AP) #function details see toolbox.r
      AV <- AP_new %*% PV
      # H_V   <- H_A_VOC(PV)[["Hn_V"]]
      # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
      # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
      H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
      H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
      # H_A_P <- H_A_VOC(AP_new)[["Hn_S_V"]]
      # H_P_A <- H_A_VOC(AP_new)[["Hn_V_S"]]
      H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
      H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
      # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
      # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
      
      # fitness function
      E_plant_new       <- F_plant(mech) 
      E_antagonist_new  <- F_antagonist(mech) 
      E_pollinators_new <- F_pollinators(mech) 
      
      # 基于优化机制判断是否保留此次基因突变
      if(E_antagonist_new > E_antagonist)
      {E_antagonis  <- E_antagonist_new
      E_plant       <- E_plant_new
      E_pollinators <- E_pollinators_new
      AP <- AP_new 
      }
    }
    
    # 更新矩阵，并计算互信息量和条件熵
    AV  <- AP %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant       <- F_plant(mech) 
    E_antagonist  <- F_antagonist(mech) 
    E_pollinators <- F_pollinators(mech) 
    
    # # 写入模拟完植食动物的 E, A_AP, A_PV
    # E[(4*n),] <- c((4*n), E_plant, E_antagonist, E_pollinators, H_V, 
    #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    E[4*n,] <- c(4*n, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
    A_AP[[(4*n)]] <- AP
    A_PV[[(4*n)]] <- PV
    A_PoP[[4*n]]  <- PoP
    
    # # 绘图
    # points((4*n),  H_P_V, col = "red", pch = 24) 
    # points((4*n),  H_A_V, col = "red", pch=  21)
    # points((4*n),  H_A_P, col = "red", pch=  22)
    # 图B
    points((4*n),  E_pollinators, col = "blue", pch = 5)
    points((4*n),  E_antagonist, col = "red", pch = 4)
    points((4*n),  E_plant, col = "green3", pch = 3)
    
    
    ## 被子植物：最大化fitness
    for (m in 1:M1) {
      PV_new <- random.sample_1element1(PV) #function details see toolbox.r
      PoV <- PoP %*% PV_new
      AV  <- AP %*% PV_new
      # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
      # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
      # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
      H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
      H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
      # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
      # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
      H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
      H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
      # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
      # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
      
      
      # fitness function
      E_plant_new       <- F_plant(mech)
      E_antagonist_new  <- F_antagonist(mech)
      E_pollinators_new <- F_pollinators(mech)
      
      # 基于优化机制判断是否保留此次基因突变
      if(E_plant_new > E_plant)
      {E_antagonis  <- E_antagonist_new
      E_plant       <- E_plant_new
      E_pollinators <- E_pollinators_new
      PV <- PV_new
      }
      
    }
    
    # 更新矩阵，并计算互信息量和条件熵
    AV  <- AP %*% PV
    PoV <- PoP %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant       <- F_plant(mech)
    E_antagonist  <- F_antagonist(mech)
    E_pollinators <- F_pollinators(mech)
    
    # # 写入模拟完植食动物的 E, A_AP, A_PV
    # E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators, H_V,
    #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
    E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators,
                   H_A_V, H_V_A, H_Po_V, H_V_Po)
    A_AP[[4*n+1]] <- AP
    A_PV[[4*n+1]] <- PV
    A_PoP[[4*n+1]] <- PoP
    
    # # 绘图
    # points((4*n),  H_P_V, col = "red", pch = 24)
    # points((4*n),  H_A_V, col = "red", pch=  21)
    # points((4*n),  H_A_P, col = "red", pch=  22)
    # 图B
    points((4*n),  E_pollinators, col = "blue", pch = 5)
    points((4*n),  E_antagonis, col = "red", pch = 4)
    points((4*n),  E_plant, col = "green3", pch = 3)
    
  }
  
  legend(x=N_sim+N_sim/40, y=0.8 , title = "Fitness", legend = c(paste0("Plant = ",E_plant), paste0("Po = ",E_pollinators), paste0("An = ",E_antagonis)),
         xpd=TRUE, bty="n",title.font = 2) 
  dev.off()
  # ggsave(
  #   filename = paste0("grid_search",i,".png"), # 保存的文件名称。通过后缀来决定生成什么格式的图片
  #   width = 7,             # 宽
  #   height = 7,            # 高
  #   units = "in",          # 单位
  #   dpi = 300              # 分辨率DPI
  # )
}








####################
# 互惠对抗1.4版
# 优化代码计算量